function [corrplanes, sups]=SLICE(data, k, fmax, minsup,maxloop)
    [M,N]=size(data);
    candiclusters={};
    candiplanes={};
    candisups=[];
    candiseeds=1:M;
    tmpdata=data;
    t=1;
    while length(tmpdata)>N && length(candiseeds)>N
        [m,junk]=size(tmpdata);
        cluster=randint(1,N,[1,length(candiseeds)]);
        for i=1:N
            cluster(i)=candiseeds(cluster(i));
        end
        [fcluster,plane,sup]=greedysearch(tmpdata,cluster,k,fmax,minsup,maxloop);
        if isequal(plane,[])==0
            candiclusters{end+1}=fcluster;
            candiplanes{end+1}=plane;
            candisups=[candisups,sup];
            %leftdataindice=setdiff(1:m,fcluster);
            candiseeds=setdiff(candiseeds,fcluster);
            %tmpdata=builddata(leftdataindice,tmpdata);
            disp(cluster);
            t
        else
            candiseeds=setdiff(candiseeds,fcluster);
        end
        t=t+1;
    end
    corrplanes=candiplanes;
    sups=candisups;
end

function [finalcluster, plane,sup]=greedysearch(data,cluster,k,fmax,minsup,maxloop)
    [datasize,N]=size(data);
    plane=[];
    sup=-1;
    finalcluster=[];
    for t=1:maxloop
        %disp(t);
        cdata=builddata(cluster, data);
        cmean = mean(cdata);
        [eigenvectors, eigenvalues]=PCA(cdata);
        [newcluster,f]=absorbpoints(eigenvectors', cmean, data, k, fmax, minsup);
        newcluster=sort(newcluster);
        finalcluster=newcluster;
        if isequal(cluster,newcluster)
            if length(cluster) >= datasize*minsup
                coe=eigenvectors(:,1)';
                const = -coe*cmean';
                plane=[coe,const];
                sup=length(cluster)/datasize;
            end
            finalcluster=cluster;
            break;
        else
            cluster=newcluster;
        end
    end
end

% Build dataset from the indices set of the data
function [cdata]=builddata(cluster, data)
    [numpoint, N]=size(data);
    numpoint=length(cluster);
    cdata=zeros(numpoint,N);
    for i=1:numpoint
        cdata(i,:) = data(cluster(i),:);
    end
end

% Absorb the points by the distance between the point and plane 
function [cluster,f]=absorbpoints(P,cmean,data,k, fmax, minsup)
    [M,N]=size(data);
    ydata=zeros(1,M);
    pydata=zeros(N,N,M);
    cluster=zeros(1,M);
    COV=zeros(N,N);
    numpoint=0;
    minicount=M*minsup;
    for i=1:M
        y=data(i,:)-cmean;
        py=P*y';
        pycov=py*py';
        v=diag(pycov);
        ydata(i)=sum(v(1:k))/sum(v);
        pydata(:,:,i)=pycov;
    end
    [junk, rindice]=sort(ydata);
    for i=1:M
        COV=(i/(i+1))*COV+(i/((i+1)*(i+1)))*pydata(:,:,rindice(i));
        V=diag(COV);
        f=sum(V(1:k)) / sum(V);
        if f > fmax
            break;
        end
        cluster(i)=rindice(i);
        numpoint=numpoint+1;
    end
    cluster=cluster(1:numpoint);
end

% Find the smallest component vector
function [eigenvectors,eigenvalues]=PCA(data)
    [M,N]=size(data);
    mn=mean(data,1);
    newdata=data - repmat(mn,M,1);
    covariance=1/M*newdata'*newdata;
    [PC,V]=eig(covariance);
    V=diag(V);
    [eigenvalues, rindices]=sort(V);
    eigenvectors=zeros(N,N);
    for i=1:length(rindices)
        eigenvectors(:,i)=PC(:,rindices(i));
    end
end

